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We present a mean-field theory for the dynamics of driven flow with exclusion in graphene-like 
structures, and numerically check its predictions. We treat first a specific combination of bond 
transmissivity rates, where mean field predicts, and numerics to a large extent confirms, that the 
sublattice structure characteristic of honeycomb networks becomes irrelevant. Dynamics, in the 
various regions of the phase diagram set by open boundary injection and ejection rates, is then 
in general identical to that of one-dimensional (ID) systems, although some discrepancies remain 
between mean-field theory and numerical results, in similar ways for both geometries. However, 
at the critical point for which the characteristic exponent is z = 3/2 in ID, the mean-field value 
2 = 2 is approached for very large systems with constant (finite) aspect ratio. We also treat 
a second combination of bond (and boundary) rates where, more typically, sublattice distinction 
persists. For the two rate combinations, in continuum or late-time limits respectively the coupled 
sets of mean field dynamical equations become tractable with various techniques and give a two- 
band spectrum, gapless in the critical phase. While for the second rate combination quantitative 
discrepancies between mean field theory and simulations increase for most properties and boundary 
rates investigated, theory still is qualitatively correct in general, and gives a fairly good quantitative 
account of features such as the late-time evolution of density profile differences from their steady 
state values. 
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I. INTRODUCTION 

In this paper we consider the dynamic evolution, as 
well as selected steady-state properties, of a generaliza¬ 
tion of the totally asymmetric simple exclusion process 
(TASEP) to two-dimensional honeycomb structures. A 
previous publication [lj focused mainly on the evalua¬ 
tion of steady-state currents for several variations of such 
structures. 

The TASEP, in its one-dimensional (ID) version, ex¬ 
hibits many non-trivial properties including flow phase 
changes, because of its collective character @i|. The 
TASEP and its generalizations have been applied to a 
broad range of non-equilibrium physical contexts, from 
the macroscopic level such as highway traffic Q to the 
microscopic, including sequence alignment in computa¬ 
tional biology 0 and current shot noise in quantum-dot 

chains ta¬ 
in the time evolution of the ID TASEP, the particle 
number nt at lattice site £ can be 0 or 1, and the for¬ 
ward hopping of particles is only to an empty adjacent 
site. In addition to the stochastic character provided by 
random selection of site occupation update |12l . 0 , the 
instantaneous current Jee+i across the bond from t to 
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£ + 1 depends also on the stochastic attempt rate, or 
bond (transmissivity) rate, pe, associated with it. Thus, 

ni( 1 — ne+i) with probability pe 
0 with probability 1 — pe ■ 

(1) 

In Ref. 11 it was argued that the ingredients of (ID) 
TASEP are expected to be physically present in the de¬ 
scription of electronic transport on a quantum-dot chain; 
namely, the directional bias would be provided by an ex¬ 
ternal voltage difference imposed at the ends of the sys¬ 
tem, and the exclusion effect by on-site Coulomb block¬ 
ade. 

Apart from the importance of generalizing fundamen¬ 
tal dynamic studies of the linear chain TASEP to higher¬ 
dimensional lattices and structures, the present work, 
and in particular its emphasis on honeycomb structures, 
is partly motivated by recent progress in the physics 
of graphene and its quasi-ID realizations, such as nan¬ 
otubes and nanoribbons 0 . Of course the TASEP, as 
described above, does not provide a realistic description 
of electronic transport in carbon allotropes under an ap¬ 
plied bias. However, in the transport context the lattice 
topology affects how currents combine, and how they 
are microscopically located, whether classical or quan¬ 
tum. It will be seen that these features show up in the 
model we treat by such effects as the sublattice struc¬ 
ture seen e.g. in steady states for the uniform hexagonal 
lattice (Secs. IllBl and IIII Cl) , and as consequent band- 
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doubling in the corresponding spectrum. An interesting 
effect, which we exploit, is the similar behavior arising in 
topologically trivial linear systems from alternating bond 
rates ISec. IIIICl) . 

Although the model is classical, so it does not display 
quantum interference effects, it is the simplest coopera¬ 
tive driven model, with consequent qualitative properties 
reflecting aspects of Coulomb blockade phenomenology in 
real experiments. 

In Ref. [l] we probed for the existence of similar spe¬ 
cific signatures by examining the behavior of steady-state 
currents for nanotubes and nanoribbons, against varying 
system sizes, and for diverse combinations of bond trans- 
missity rates, as well as distinct sets of boundary condi¬ 
tions along the flow direction, namely periodic (such as 
to make the system ring-like) and open (with assorted 
values for injection and ejection rates at the ends, to be 
recalled in detail below). 

The latter case of open systems, with open boundary 
conditions at the ends, is by far the most challenging, 
richest, and most illuminating one, so it (alone) is the 
case here considered. As in Ref. [lj, the present study 
makes complementary use of mean field analysis and nu¬ 
merical simulations. 

In Section El a mean field theory is presented which 
describes the time evolution of ensemble-averaged site 
occupations under TASEP rules, and applies both to the 
two-dimensional structures under specific consideration 
here and to their linear chain counterparts. Section nn 
deals with numerical tests of the theory given in Sec¬ 
tion El In Section IIVI we summarize and discuss our 
results. 


II. MEAN-FIELD THEORY 

For analytic tractability we shall only consider cases 
where mean flow direction is parallel to one of the lattice 
directions, and bond rates are independent of coordinate 
transverse to the flow direction. These configurations 
have no bonds orthogonal to the mean flow direction; 
thus they fall easily within the generalized TASEP de¬ 
scription to be used, where each bond is to have a definite 
directionality, compatible with that of average flow. 

Also, we consider structures with an integer number 
of elementary cells (one bond preceding a full hexagon) 
along the mean flow direction. See Fig. [T] 

From Ref. [l] we have to expect a two-sublattice charac¬ 
ter in general, each being of similar character to those for 
chains. For a special choice of the bond rates defined in 
Fig. [T| [p = 2 q, see the discussion of Eqs. © © below] 
the steady state sublattices reduce in mean field to that 
of an equivalent uniform-rate chain [I[. 

Throughout this paper only axially symmetric bound¬ 
ary conditions will be considered, and no rate disorder 
will be allowed for. Then, in general, the (mean) dy¬ 
namic configurations are translationally invariant in the 
direction transverse to the tube axis. Consistently with 
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Figure 1. Schematic sections of a nanotube, showing (top 
to bottom): injection region, midsection, and ejection region. 
Average flow direction is from top to bottom of the figure. 
Bond rates are p = 1 for bonds parallel to average flow di¬ 
rection, q otherwise. Injection (a) and ejection (/?) rates are 
shown next to corresponding (injection and ejection) sites. 
Periodic boundary conditions across are omitted for clarity. 


this, we denote the average occupations at sites labelled 
by the longitudinal coordinate l (1 < £ < N ) by x(£,t) 
and y{£, t ) with t odd and even respectively, correspond¬ 
ing to the two sublattices (see Fig. lU). 

Using mean field factorization, the currents on the two 
different types of bond are 


Jee+i =pxe{ 1 - ye+i) 

{£ odd) 

(2) 

I<ee+i = qye(l- xe+i) 

(£ even) . 

(3) 

Then the general equations for xe, 
are 

ye at interior 

sites £ 

xe = 2Ke~ie — Jee+i 

(£ odd) 

(4) 

ye = Je-u — 2Kee+i 

(£ even) . 

(5) 

From boundary injection and ejection at sites £ = 

: 1 and 


N, both on the x —sublattice (£ odd), incoming and out¬ 
going currents are 

a(l-£i) = Ji (6) 

2 K n = x N p . (7) 

In the steady state where xe = ye = 0, all £, these discrete 
equations specify discrete current balance, making Jee+ 1 
and 2Kee+\ equal and bond-independent (= J, say), and 
making xe and ye reduce to steady state values xe, ye , 
where 


a (1 — X\) = J = /3 Xn ■ 


( 8 ) 
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The distinct steady state sublattice characteristics are 
seen in the (in general) distinct t— dependent density pro¬ 
files xg, yi which are provided by Mobius map relation¬ 
ships between if and ye+i resulting from specified J and 
K (= J/2). 

From Eqs. © ©, it is easy to see (and was exploited 
in Ref. |l|) that the sublattice distinction goes away for the 
special case p = 2 q. Here the nanotube steady state is 
that of an equivalent linear chain, having density profile 
in general with tanh or tan dependences on t. 

The value of a continuum approach to the mean field 
dynamics of the uniform linear chain is well known 0,1 
Il3j |. and it exploits a linearization of the continuum mean 
field dynamic equations using the Cole-Hopf transforma¬ 
tion [3, 0. We show in Sec. Ill Al that this technique 
can also be successfully used for the nanotube with rates 
2q = p = 1 (for convenience) and axial symmetry. 

In Ref. [l] it was shown that for the general case p 2q, 
Mobius maps still apply, from which steady state density 
profiles are again predicted to be of tanh or tan form, but 
in general different on the two sublattices. Even though 
on each sublattice separately continuum viewpoints can 
still apply (e.g. not too far from critical conditions), stan¬ 
dard Cole-Hopf transformations no longer linearize the 
coupled nonlinear dynamic equations. Nevertheless, in 
Sec. HH1 (i) we are there able to use another lineariza¬ 
tion procedure, on the discrete equations for the dynam¬ 
ics, which gives an asymptotically exact representation of 
the mean field dynamics at very late times; and (ii) fur¬ 
thermore, we can exploit arguments (see Appendix E} 
based on the existence of two separate relaxation time 
scales, from which it follows that a continuum-like pic¬ 
ture is in fact feasible for not very short times. It will 
be seen that this, combined with simulation, can give a 
particularly clear and direct probe of critical dynamics. 


A. Continuum approach for p = 2q = 1 


In this case, no longer needing to distinguish sublat¬ 
tices, the notation p(£, t ) can now be used for the density 
profile. The continuum version of the bond current is 
then 


from which one arrives at the following form of the steady 
state profile: 

P = \ + \ Z tanh ~ £ °)] > ( 10 ) 

from J = constant = (1 — Z 2 )/4, with Z real or pure 
imaginary depending on whether the steady state current 
is less or greater than the critical value J c = 1/4. The 
resulting continuum dynamic equation 


dp 

at 


d_ 

d£ 


P { 1 - P) 


1 dp 

2 d£ 


( 11 ) 


is easily reduced to a linear (diffusion) equation for the 
variable it, by the Cole-Hopf transformation [15,, 16] 


1 1 a 1 
"-2 = 2m *"“■ 

(12) 

A general solution reducing as t — > oo to the steady state 
profile p given in Eq. (HOD is 

it = it + E , where 

(13) 

u = cosh [Z(£ — £o)] e = z 4 

(14) 

E = ^ (af e ^ + a_£ e - ^) 4 

(15) 


C 


where the sum is over (/s, in general complex, satisfying 
5ft C 2 < 5ft Z 2 . For the validity of the continuum approx¬ 
imation, Z and all effective £'s arising should be small. 
The boundary conditions, Eqs. © and (T7)l . which deter¬ 
mine them can be rewritten as (for all t ) 

p(0, t) = a (16) 

P (N + l,t) = 1-/3, (17) 

where p(0,f) and p(N + 1 ,t) are the extrapolations of 
the solution, Eqs. (ED ED, of the dynamic equations 
to the fictitious sites immediately outside of the system 
boundaries. These then have to be satisfied by the steady 
state part p = (1/2) + (l/2)d\nu/d£, as well as by the 
time-dependent parts of the extended p. The require¬ 
ments on p give Z , £q, in particular requiring Z real for 
a < 1/2 or [3 < 1/2; or Z pure imaginary for a > 1/2 and 
/3 > 1/2. From the time-dependent parts the boundary 
conditions then require d In E /d£ equal to p\ = 2a — 1 
and /i 2 = l — 2/3at£ = 0 and £ = L = N+ 1 respectively. 
That leads to 


a -C = C ~ Pi = e 2 Cl f C ~ P'2 \ 
®c C + Mi vC + M 2 / 


(18) 


giving both the allowed complex wave vectors £, and the 
ratio of associated amplitudes. Initial conditions then in 
principle complete the determination of all amplitudes, 
by the analogue of Fourier analysis. 

Some special cases will be of interest in what follows, 
namely a = (3 and a + (3 = 1. 

For a = /3, the open boundary condition restrictions 
make £q = A/2, and, for a = /3 < 1/2, Z is real, say Z = 
K, with K = 1 — 2a + 0(e ~( 1-2a ) L ) - so the dynamics 
is relaxation to the steady state of the low current phase, 
having a kink in the middle of the system; while, for 
a = (3 > 1/2, Z is pure imaginary, say Z = iQ, with 
Q = (2/A)((7t/ 2) — 7r/[A(2a — 1)]), and the relaxation is 
towards the high current phase steady state. 

For the critical subcase a = f3 = 1/2, one has Z = 0, 
( n = mri/L = iq n , a_£ = = a n . So for this case it = 0 

and u = E = a n (e^ n + e -< »") * making 

(p t \ = l_ 1 E«=i sin (q n £) 

22 


(19) 
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A given initial profile p(£,0) would complete the de¬ 
termination of p(£, t ) by providing the coefficients 
a n , by the equivalent of Fourier cosine analysis of 


exp 


f 0 d£' (2p(£', 0) — 1)) 


in the present case. 


For an initially empty lattice, for example, this gives 


_ 2 [1 - (-!)» exp(—L/2)] 

L 1 + Q-n 


2 

L 


[i + qI 


i -i 


( 20 ) 


Then, for late times t > ( L/tz ) 2 , 


1 

1 . 

£n£\ 


i ( 

rr \ 2 

pra 2~ 

2 9i sin 


exp 

2 (" 

l) 1 _ 


while for early times 1 <Cf < (L/n) 2 


( 21 ) 




Wt)= / dC 


cos((£)e sC 2 * 

1 + C 2 


( 22 ) 

making I(£,t)y/t essentially a function of £/y/t, and p 
linear in £ (p « (1/2) — (£/21)) up to £ ~ 0(y/i). This 
is of course related to the buildup of density from the 
injection site, and is evident in simulation results shown 
in Sec. IIII1 see Fig. [5] 

For a + f3 = 1, the boundary restrictions on the steady 
state are consistent with Z = 1 — 2a = A, and £o oo 
for A > 0 and £q —>■ —oo for A < 0 (kinks far outside 
of the system). From the other boundary restrictions, 
( n = niri/L , and a-c n / a Cn = (£„ + \)/(( n — A). These 
make the steady state u proportional to exp [— X£ + ^A 2 t] 
and 


E 

u 


S = ^ (a£ n e^ + a_£ n e < ’ r,e ) 

71=1 


e M e i((Z~* 2 )t 


(23) 


Then the time dependent density profile becomes 


1 dS/d£ 

21 + 5 ' 
(24) 

In this case the relaxation is towards the constant (fac- 
torizable) steady state profile p£ = a; at very late times 
one has 


1 1 d . 1 1 d , 

p= 2 + 2S 1,, " = 5 + 2af 1, ' [ " (1+s,1 = “ + 


p — a 


IdS 
2 d£ 


£(*-t 


71=1 

x exp 




(A 2 + ( mr/L ) 2 ) t 




x 


(25) 


where A = 1 — 2a. For a = 1/2, Eqs. ¥IM and (051) reduce 
to Eq. d). 

The late-time results in Eqs. m and (E3 above, and 
others to be given in Sec. IllBl [especially Eqs. (l47l) and 
(BED] can give guidance beyond the mean field regime 
used to obtain them. The correspondence, within mean 
field, between chain and nanotube for the case p = 2q 
(given for the steady state in Ref. |l] and extended here 
to dynamics) implies the same mean field exponents, and 


this is seen also for 2 q ^ p below, see Sec. Ill Bl In partic¬ 
ular the functional dependences on £/y/t and t/L 2 seen 
above [ in Eq. CD, and in the equivalent Eq. C3 for 
A = 0] correspond to the mean field value 2 of the dy¬ 
namic critical exponent z. But one can reasonably expect 
the (wide) nanotube to have different critical exponents 
from those known for the chain, since the two have dif¬ 
ferent dimensions. 

The simulation method in Sec. nn is able to exhibit 
these differences, and the mean field analytic results sug¬ 
gest a direct method to find them, by exploiting the late 
time behavior, in particular by using the slowest-relaxing 
mode. 

The results in Eq. CD and (1251) (the latter, from 
just the n = 1 term) provide mean field examples of 
that mode, and suggest that its isolation, by working 
at late times, particularly when the system is relaxing 
to a uniform steady state [ using p(£, t) — pi\, can give 
the most unencumbered way of numerically investigat¬ 
ing the critical dynamics. Finite-size scaling using fitting 
forms for p{£,t) — p{£), like in Eq. (HHll or in the n = 1 
mode of Eq. CD, but with the time-dependent factor re¬ 
placed by exp[— ctL~ z ] are suggested: the general form 
f(£/L,t/L z ) could, from the last surviving eigenmode 
of the evolution operator e~ Ht , go over to a factoriz- 
able form having an e~ t !' r time-dependent factor, with 
r ~ L z , and a spatially-dependent factor with nodes 
near £ = 0, L (from boundary conditions) and a sym¬ 
metric form [ like in Eq. (ET1) ] or with an extra factor 
as in Eq. CD, the latter in cases with p ^ 1/2. These 
ideas are exploited in Sec. ED both for the chain and for 
the nanotube. 


B. Discrete late-time method, for p 2q 

Here we develop an analytic method for the late time 
dynamics, which is applicable for general rates a, /3, p, 
q where sublattices are distinct and remain so even in 
the eventual steady state. Unlike Sec. Ill Al using the 
continuum approach, the method proceeds from the dis¬ 
crete mean field dynamic equations and linearizes them 
by working to first order in differences of site densities 
from steady state values. 

The discrete steady state densities are determined 
by the Mobius maps introduced in Ref. [l|, which re¬ 
sult from steady state internal current balance, together 
with boundary conditions, as explained after Eq. ([7[l. If 
these densities are site-dependent the difference dynam¬ 
ical equations resulting from the linearization procedure 
have site-dependent coefficients, making them in general 
intractable. For 


a = 2g(l -p) (p= 1) (26) 

the steady-state densities given by the Mobius mappings 
can be uniform on each sublattice, while in general re¬ 
maining distinct. 

















5 


The analysis now to be given treats that case, at gen¬ 
eral q , for which the coupled linear difference equations 
have constant coefficients. Their solutions are linear com¬ 
binations of factorizable solutions, involving a secular 
relation between the frequency and complex wave vec¬ 
tors involved. The boundary conditions determine the 
allowed values of the complex wave vectors and relation¬ 
ships between amplitudes of degenerate components. 

The uniform steady state density profile values x, y 
on the two sublattices correspond to fixed points of the 
discrete Mobius maps. Such fixed points are directly 
available from the basic internal and boundary current 
balance equations 

a(l — x) = x{\ — y) = 2qy{l — x) = [3x . (27) 

Provided a = 2q(l — (3) these result in 

x = —jtr ; y = 1 - £ • ( 28 ) 

a + p 

An important subcase to be distinguished and developed 
later in this section is the critical situation, where the 
two fixed points for each sublattice Mobius map coincide 
(corresponding to Z = 0 in the continuum steady state 
description in Eq. m, see Sec. EM- 

Starting from the discrete mean field dynamical 
Eqs. (U) and ([5]) the linearization procedure, valid for 
sufficiently late times, is implemented by inserting xg = 
x + Sg, yg = y + £g and including only terms up to first 
order in Sg, Eg. 

The zeroth order terms involving only x and y are those 
appearing in the steady state current balance, so they 
cancel. The resulting coupled linear difference equations 
for the time-dependent <5^, eg are solved by superpositions 
of factorizable solutions of the form 

Sg = gc, exp(C£ - A t) (29) 

eg = exp(££ — A t) (30) 

for specific £- and A-dependent ratios h^/g^ provided £ 
and A satisfy the secular relation 


Then, with ry = — (p, the degeneracy condition becomes 

r )i = —r] 2 - That allows the superposition of degenerate 
modes for Sg to be written as 

Si = E ^ + 9*-v e~ ve ) e* e~^ + ^ , (35) 

v 

and similarly for eg (where h^q replace g<f,±q). 

The secular relation between A and £ can be rewritten 
as one between A and g using 

S(C = V + <p) = So — y{g) where 

y(g) = (e” + e-*) . (36) 

For the boundary conditions to be maintained by the 
full time-dependent profiles xg = x + 5g, yg = y + £g, 
the differences Sg, Eg have both to vanish at i = 0 and 
£ = L at all times. That requires g^+q + g<f>-q = 0 = 
h(f, + q + hcji-q and e 2r]L = 1, so the allowed g's are g n = 
imi/L = iq n . 

Consequently the space- and time-dependent sublat¬ 
tice density profiles are, to linear order, 

xg(t) = x + ^ G n smq n £e^ e e~ Xnt (37) 

n 

yi(t) = y + F H n sin q n t e^ e~ Xnt (38) 

n 

where 



and A„ satisfies 

A 2 - r\ n + S 0 - y(iq n ) = 0 (40) 

where r and Sq are given in Eqs. and (E£T1) . and 

y{iq n ) = 4g y/xy (1 - x) (1 - y) cos q n , (41) 


A 2 - rA + S{ C) = 0 , (31) 

where 

r = 1 + 2q + (1 — 2 q)(x - y) ; 

^(C) = 5o-(5 + e c + 5_e- c ) (32) 

with 

S 0 = 2q(l -x-y)+ 4 qxy 
S + = 2 qxy 

S- = 2q(l-x)(l-y) . (33) 

To fit the boundary conditions at all times it is necessary 
to combine degenerate modes, i.e., modes with ^ £2 
such that A(£i) = A^)- A sufficient condition for this is 
5(Ci) = 5(<^2), from which 


e C 1+ C 2 = |^^ e 2* . 
* 3 + 


(34) 


with x. y given by Eq. (1251) . 

The coefficients G n and H n (piig^-q and 2 ih^+q, re¬ 
spectively) are in principle determined by initial states. 
For initial states ££(0), yg(0) in the linearization regime, 
they are the coefficients in the Fourier sine series for 
xg(0) — x, yg( 0) - y respectively. 

The very late time behavior, from the decay of the last 
surviving time-dependent mode, is described by xg(t) — x, 
yg(t) — y both proportional to sin(7r£/L) e _Alt , with 
<p from Eq. (1551) and Ai = \ [r— yjr 1 — 4(So — y{iqi))\ 
from Eqs. (EHD dSD • 

In general, the distinct sublattices give rise to a two- 
branch spectrum, which makes the late-time dynamics 
for the cases with 2 q ^ p very different from that with 
2 q = p discussed in Sec. Ill AI The spectrum is in general 
gapped even in the infinite-system limit (hin^^oo Ai >0) 
as a consequence of non-zero <)>; the gap goes away (as 
<f> —> 0) only in the critical cases, which we now discuss. 
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The critical steady state has constant (coincident fixed 
point) values x*, y* for x, y, related to a critical cur¬ 
rent J c on the bonds with rate p, and to critical bound¬ 
ary rates (a c ,/3 c ) by current balance equations of type 
Eq. (1271) . where each current is J c such that the cor¬ 
responding sublattice Mobius maps each have coincident 
fixed points. Withp = 1, that requires [J c (l— 2q) —2q] 2 = 
16 q 2 J c , which makes x* + y* = 1, hence 

<t> = 0 , S 0 =2S+=2S- = 2qx*y* . (42) 

That in turn makes 


S(C)=S{rj)=4qx*y*( 1-coshC) (43) 


and the development in Eqs. ( 1371 ) ( 1711 ) simplifies. The 
results for the time-dependent density profiles become, 
to linear order, 

xe(t) = x* +'^2G n sinq n £e~ Xnt (44) 

n 

yi{t) = y* + y^ j H n sin q n i e~ Xnt (45) 

n 

where 


A„ = - r±\Jr 2 - 16qx*y*(l - cos q n ) = A ±(q n ) ■ 

( 46 ) 

So, (f> = 0 has produced a gapless spectrum in infinite 
system limit, for the critical system, and we now have 
the analogue of acoustic and optic modes. 

For the finite critical system, the very late behavior 
of the profiles on each sublattice is (using the slowest 


xg(t) = x* + Gsin 
ye{t) =y* + H sin 


1, with A_) 


/L)t 

L 

(47) 

— p-A-O /L)t 

L 

(48) 


with 


. (7 r\ S(C) 4qx*y* / 7 r \ / tt\ 2 

7 2 ~-V l ( 1 -“ s zMz) ■ 

(49) 

The condition a = 2q(l — 0) for uniform steady state 
densities, presently applying, reduces for q = 1/2 to a + 
/3 = 1, which is a case discussed for general A = 1 — 2a 
in Sec. Ill Al That case becomes critical for A = 0, and 
then the present formulation with <fi = 0 together with 
the resulting Eqs. (1721) (1771) all apply to it (reproducing 
results in that Section). 

A particular important special case is that for the 
uniform-rate nanotube, where p = q = 1 and x* = 2—\f2 1 
y* = s/2- 1, a c = 2(s/2- 1 ), 0 = 2-02 (from Ref. HI), 
agreeing with Eqs. EH and (1281) . 

The distinction, one or two bands (from 2 q equal to p 
or not) is a special feature of the nanotube coming from 
its possible sublattice character, and shared with the 
TASEP chain with alternating bond rates p , 2g, which 
has equivalent mean field steady state and dynamics. 


III. NUMERICS 


With open boundary conditions at the ends, a nan¬ 
otube with N r elementary cells parallel to the flow direc¬ 
tion, and N w transversally, has = N w x (4N r + 1) 
sites and = N w x (6N r + 2) bonds (including the 
injection and ejection ones). 

When dealing with strictly ID geometries, for ease of 
pertinent comparisons with nanotubes we generally took 
systems with a number of sites N = 4M +1, M being an 
integer. 

Here we shall only use so-called bond update proce¬ 
dures, defined in Ref. Q] and briefly recalled below. For 
a description of the closely-related site update process, 
and pertinent comparisons with bond update, see Ref. [J 

For a structure with Nb bonds, an elementary time step 
consists of Nb sequential bond update attempts, each of 
these according to the following rules: (1) select a bond 
at random, say, bond ij, connecting sites i and j\ (2) if 
the chosen bond has an occupied site to its left and an 
empty site to its right, then (3) move the particle across 
it with probability (bond rate) pij. If the injection or 
ejection bond is chosen, step (2) is suitably modified to 
account for the particle reservoir (the corresponding bond 
rate being, respectively, a or 0). 

Thus, in the course of one time step, some bonds may 
be selected more than once for examination and some 
may not be examined at all. This constitutes the random- 
sequential update procedure described in Ref.fl2l which is 
the realization of the usual master equation in continuous 
time |l2]. In our simulations, the goal for ID uniform sys¬ 
tems is to have numerically-generated profiles approach 
the exact steady-state ones given by the operator alge¬ 
bra described in Ref. [5, which are an important baseline 
in our work and, as recalled in Ref. Il2l . correspond to 
random-sequential update. For consistency, and ease of 
comparison between different sets of results within the 
paper, we also use random-sequential update for all other 
cases, namely honeycomb geometries and non-uniform 
ID systems. Note that other types of update are pos¬ 
sible (e.g., ordered-sequential or parallel), the resulting 
steady-state phase diagrams in ID being similar in all 
cases (but not identical: even the average stationary cur¬ 
rent differs in either case, see Table 1 in Ref. Il2ll . 

For specified initial conditions, we generally took en¬ 
semble averages of local densities and/or currents over 
10 6 10' independent realizations of stochastic update 
up to a suitable time f m ax, for each of those collecting 
system-wide samples at selected times. 

For uniform ID systems and nanotubes with p = 2 q, 
the exact steady-state density profiles {xe}, known in 
ID for any a, 0 and N [5( are used as a baseline 
from which to subtract our late-time simulational results 
{xe(t)}, thus focusing on the evolution of difference pro¬ 
files 5x0t) = x0t) — X£. For nanotubes with p 2 q, or 
chains with non-uniform rates, both cases considered in 
Sec. 1III Cl no such guidance is available. One must then 
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Figure 2. Linear chain with N = 41 sites (L = N + 1), 
a = P = 1/2. Double-logarithmic plot of (negative) ini¬ 
tial slopes ( S ) of short-time density profiles against time t 
(points). Continuous line is the mean-field prediction S = 
1/(2 1), see Eq. (1501) . The vertical dashed line indicates the 
approximate limit of validity of the short-time regime (see 
text). 

resort to numerically-generated steady state profiles. 

A. p — 2q, a = /3 = 1/2 

We started by checking the predictions given in 
Sec. Ill Al for the time-dependent density profiles of a ID 
system starting from an empty lattice. Eq. (1221) predicts 
that for short times t <C ( L/tt ) 2 , 

(50) 

near the injection edge, up to £ ~ 0(y/t). For a 
chain with N = 41 sites, we evaluated the initial slope 
dp/d£\\g=o at assorted short times, from straight-line fits 
of ensemble-averaged densities at the three leftmost sites. 
Results are shown in Fig. [2] One sees that agreement 
between theory and numerics is rather satisfactory, espe¬ 
cially if, drawing on the last two paragraphs of Sec. Ill Al 
and on previous knowledge of the anomalous scaling for 
ID systems at (a,/3) = (1/2,1/2), one restricts oneself 
to data for t < ( L/t :) 3 / 2 [ as opposed to t < (L/n ) 2 from 
the mean-field picture leading to Eq. (l22ll ]. 

Next we checked the late-time behavior, both for ID 
systems and for nanotubes. Fig. [3] shows a fit of Eq. (BT1) 
to the ensemble-averaged density profile for a ID system, 
starting from an empty lattice at t = 0. While the quality 
of fit is good, with y 2 per degree of freedom (Xd of ) equal 
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Figure 3. Linear chain with N = 29 sites, a = /3 = 1/2. Plot 
of late-time density profile, starting with an empty lattice at 
t — 0. Continuous line is the fit to a sine form, see Eq. in 
and text. 

to 1.35, one sees that small systematic deviations still 
remain near the left (injection) edge. Going over to later 
times in order to evince the suppression of such deviations 
would necessitate much narrower error bars (since one 
would be analyzing profiles much closer to the asymptotic 
regime), and consequently much longer simulations, than 
in our current setup. 

Nevertheless, we now show that it is possible to extract 
rather accurate estimates of the dynamic exponent z from 
our data in present form, by once again referring to the 
ideas sketched in the last two paragraphs of Sec. Ill Al 
Specifically we rewrite Eq. m as 

p&t) = \ exp { _c ^} ’ ( 51 ) 

i.e., while assuming factorization of the 0. and t depen¬ 
dences, we allow z to be a variable parameter. For 
fixed L and a set of suitable t values, fitting numerically- 
generated profiles to the sine dependence in Eq. m pro¬ 
duces a sequence of estimates of 

a*(L,t) = a'(L) expj-c-^j ; (52) 

the latter set is then fitted to 

a*(L,t) = ao(L) exp{— c'(L) t} , (53) 

with ao(£), c'(L) as fitting parameters. Finally, varying 
L one fits the corresponding sequence of c'(L) to a power- 
law in L, thus extracting z. 









































We proceeded as just outlined for: (i) ID systems, 
starting with an empty lattice; (ii) ID systems, starting 
with a "sine-like" profile, i.e., 


nt( 0) 


1 £ < f or l > 

0 f <£< ^ 


(54) 


in order to check how sensitive the small late-time sys¬ 
tematic deviations, referred to above, were to the choice 
of initial condition; (iii) nanotubes with N w = 14 elemen¬ 
tary cells across and varying length Ay; and finally (iv) 
nanotubes with N w = N r cells, i.e. aspect ratio equal 
to unity. In the latter two cases, sine-like initial profiles 
were used. 

For (i)-(iii) we took N = 29, 41, 53, and 69 (corre¬ 
sponding, for nanotubes, to Ay = 7, 10, 13, and 17) and, 
for each of these, five N- (or independent values of t in 
the late-time approach to steady state. We found that 
using a sine-like profile as initial condition does slightly 
improve the quality of profile fits to Eq. (I2D- For ex¬ 
ample, in the corresponding case to that illustrated in 
Fig. H we found x t 2 )of = 0.91, about a third less than for 
an enrpty-lattice start. 

By following the fitting procedures delineated above 
our final results were 2 = 1.51(1) in case (i), z = 1.54(1) 
in case (ii). The main diagram in Fig. U illustrates how 
well the numerically-evaluated coefficients a*(L , t) follow 
an exponential decay in time. That, as well as the smooth 
power-law fit of c'{L) against L shown in the inset, gives 
strong support to the ansatz described in Eqs. m (El- 

Analysis of case (iii) for the nanotube produced a less 
clear-cut picture concerning the final estimate of z. Al¬ 
though the exponential decay in time of the a*(L,t) 
still holds to excellent accuracy, resulting in the coeffi¬ 
cients c'(L) listed under the heading (a) N w = 14 in Ta¬ 
ble 1 a single power-law fit of the latter against L gives 
z = 1.76(2). By drawing on ideas for successively iterat¬ 
ing sequences of finite-size approximants of quantities of 
interest we produced a set of two-point fits of data 
for pairs {L U L 2 ) = (30,42), (42,54), and (54,70). Plot¬ 
ting such set against 2/ (L 1 +L 2 ), we arrived at the follow¬ 
ing extrapolated values for 2/(Li + L 2 ) —> 0: z = 1.58(1) 
for a linear fit, 0 = 1.51(2) for a parabolic fit, see Fig. [5] 

In case (iv) we took N r = N w = 8, 12, 16, and 22. The 
sequence of coefficients c'(L), obtained along the same 
lines already described, is given in Tabled] under (b) As¬ 
pect Ratio= 1. As shown in Fig. [5] by iterating two-point 
fits for pairs of successive lengths ones gets an increas¬ 
ing sequence of estimates of z against increasing L. A 
straight-line fit gives an extrapolated z = 2.04(4). So this 
indicates that, while keeping N w > 1 fixed one gets essen¬ 
tially one-dimensional (critical) behavior, allowing for a 
constant aspect ratio of order unity one picks (asymptot¬ 
ically) the true two-dimensional dynamics. Furthermore, 
numerics indicate that the latter is characterized by the 
mean field exponent z = 2. 

Going back to the data for fixed N w , for the nanotube 
with p = 2q = 1, (a, (3) = (1/2,1/2) there appears to 



Figure 4. Main diagram: log-linear plot of a*(L, t) of Eq. (1521) 
against t for linear chain with N = 29 sites. The continuous 
line connects numerically-obtained points. Initial condition: 
sine-like. Inset: double-logarithmic plot of c'(L) of Eq. (15311 
against L = N + 1. The continuous line is a fit of data to 
c'(L) ~ L ~ z , with 2 = 1.51. Initial condition: empty lattice. 


Table I. For nanotubes with p = 2q = 1, (a, /3) = (1/2,1/2), 
late-time coefficients c'(L) of Eq. (1531) . obtained by the fitting 
procedure described in the text, for varying system lengths L. 
(a): fixed width N w = 14 hexagons; (6) aspect ratio = 1. 


L 


c\L) 

(a) N w = 14 

30 


0.00859(12) 

42 


0.00458(2) 

54 


0.00291(3) 

70 


0.00185(2) 

(6) Aspect Ratio = 1 

34 


0.00670(2) 

50 


0.00340(3) 

66 


0.00203(1) 

90 


0.00113(1) 


be a slow crossover towards 2 = 3/2 behavior against 
increasing system size, which does not have a parallel in 
strictly ID systems. 

We have checked this scenario by investigating a 
steady-state quantity which is well-known to display sig¬ 
natures of anomalous scalin g, n amely the cumulants of 
the integrated current fill . I18j . Denoting by J the 
steady-state average current through a specified bond, 
say the one linking sites £ and £ + 1, and Jee+i(t') its 
instantaneous value, the associated integrated charge is 
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2/(L[ + L 2 ) 


Figure 5. Nanotube with p = 2q = 1 at critical point 
(a, {$) = (1/2, 1/2). Points are estimates of dynamical expo¬ 
nent 3 resulting from two-point fits of c'(L) in Table|I]for pairs 
of successive lengths (Li,Z/ 2 ), against 2/(Li + 1/2). Squares: 
fixed width N w = 14. Circles: fixed aspect ratio (A. R.)= 1. 
Full lines: linear fits. Dashed line: parabolic fit [ for N w = 14 
only] (see text). Initial condition: sine-like in all cases. 

Qei+i(t) = fo Jee+i(t') dt'. Usually one removes the lin¬ 
ear term, and considers 

Q{t) = Q(t ) - Jt , (55) 

so (Q(t)) = 0. For ID TASEP at (a,/3) = (1/2,1/2) the 
second-order cumulant (( Q 2 )) of the integrated current 
has been shown fill.[ill] to exhibit anomalous scaling, i.e., 
(( Q 2 (t))) ~ t 1//z with 2 = 3/2 along a time "window" of 
width determined by system size ("normal" scaling would 
correspond to (( Q n {t))) ~ t for all n). In Fig. [G] we 
show data for both ID systems, and for a nanotube with 
N v , = 12, N r = 10 (N = 41). The apparent behavior 
oc t 0 - 57 exhibited for 200 < t < 5 x 10 4 by the latter 
is consistent with {{Q 2 (t))) ~ f 1 / 2 , using the effective 
exponent z = 1.76(2) found from a global analysis of the 
c'(L) for fixed N w of Table [T) 

Still for the nanotube, one can see behavior compatible 
with (( Q 2 (t ))) - t 2 / 3 for 5 x 10 4 < t < 2 x 10 5 , until it 
crosses over to "normal" scaling (( Q 2 {t))) ~ t (of course 
the latter also takes place for ID systems, see the cor¬ 
responding data in Fig. [G]). The narrowness of the t 2 / 3 
"window" is most likely related to the relatively small 
(longitudinal) system size N Eua. 

So in the quasi-one dimensional limit for the p = 2q 
nanotube at criticality, the evidence provided both by dy¬ 
namics (from the scaling of the c'(L) of Eq. (IGG1) against 
L) and steady-state (from the scaling of (( Q 2 {t))) against 
t ) consistently points to an apparent z — 1.76 for rela- 
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Figure 6. Points represent numerically-evaluated second cu¬ 
mulant ({Q 2 (t)}) of integrated steady-state current versus 
time t, for (a, /3) = (1/2,1/2). ID: linear chain, N = 600 
(adapted from Ref. Il8). NT: nanotube of width N w = 12 
hexagons, N = 41, with bond rates p = 2q = 1. Lines in¬ 
dicate power-law dependence with exponents as shown (see 
text). 

tively short systems, and/or short times (the latter, af¬ 
ter full onset of the steady-state regime), followed by a 
crossover towards 2 = 3/2 in this case. 


B. p = 2 q, a + /3 = 1 

For a + (3 = 1, away from the critical point which 
was the subject of Sec. IIII Al we took a point in the low 
current phase of ID TASEP, namely a = 0.3, /3 = 0.7. 

Considering ID systems, starting from an empty lat¬ 
tice, we adapted Eq. (l?5l) for very late times such that 
only the n = 1 term in that Equation still survives. In or¬ 
der to investigate density profiles in this regime we write: 

p(£, t ) = 0.3 — a(L, t) sin i (56) 

where a(L, t) incorporates the exponential time depen¬ 
dence in Eq. (1251) . and the factor b in Eq. (IGG1) is pre¬ 
dicted to be b = A = 0.4. In Fig. [7J for ID TASEP 
with N = 29, curve (I) [full red line] shows the best 
fit of Eq. (1561) to the simulational results given there, 
corresponding to a = 35(4) x 10 -5 , b = 0.144(6), with 
Xj of = 1.8. Curve (II) [dashed blue line] is the prediction 
of Eq. (EGl) with A = 0.4, with the {aj„} adjusted to an 
empty-lattice initial condition and using only the n = 1 
term. 
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Figure 7. Linear chain with N = 29 sites, a = 0.3, /3 = 0.7, 
p = 2q = 1. Plot of late-time density profile, starting with 
an empty lattice at t = 0. Curve (I) is a fit to a sine plus 
exponential form, with a and b of Eq. ([56| as adjustable pa¬ 
rameters; curve (II) is the prediction from Eq. (1251) . adjusted 
to the empty-lattice initial condition, and using only n = 1; 
see text. 


Although their overall shape is similar, curves (I) and 
(II) significantly differ in (a) the depth and, to a lesser 
extent, location, of the minimum on the right-hand side, 
and (b) the nearly-horizontal segment stretching almost 
midway through the system, exhibited by curve (II), 
which has no counterpart in curve (I). While making 
t ft 110 in Eq. (I251l . instead of "simulation time" t = 120 
reproduces the minimum value shown by numerical data 
(its location, however, remaining unchanged within one 
lattice spacing), point (b) is a permanent feature of the 
theoretical prediction which reflects the large value of 
A = 0.4 in the profile’s exponential ^—dependence in 
Eq. (l25l) . 

The discrepancy between the optimally adjusted value 
of the exponential prefactor b of Eq. (ESP, on the one 
hand, and the theoretical prediction of A = 1 — 2a on 
the other, is undoubtedly significant. This indicates that, 
although simple adaptations enable it to give an accurate 
description of the critical systems of Sec, fill Al the mean- 
field theory given above does not quantitatively account 
for the effects of a characteristic inverse length A ^ 0 in 
a similarly straightforward way. We have found |19| that 
a formulation including the effects of stochastic domain- 
wall hopping [20l [2^ | on early- and late-time profiles can 
account for most of the quantitative mismatches between 
mean-field theory predictions and simulational results for 
non-critical cases. 

However, in the present work we limit ourselves to 


analysing the extent to which the mean field theory of 
Sec. HU can provide useful clues to the actual behavior of 
numerically-generated samples. Thus, here we attempt 
a procedure similar to that described in Sec. IfflAl for 
extraction of the dynamical exponent. 

In addition to ID systems, and similarly to Sec. IIII Al 
we considered nanotubes both (1) with N w = 14 ele¬ 
mentary cells across and varying length N r . and (2) with 
unit aspect ratio (N w = N r cells). The time dependences 
predicted respectively in Eqs. ED, related to the critical 
("gapless") phase and (E5l) for the "massive" or "gapped" 
phase, differ in that the decay rate in the latter has an 
L —independent term, the gap [equal to A 2 /2 ], related to 
the characteristic inverse length A. 

In an attempt to give similar relative importance, when 
compared to the gap contribution, to the finite-size de¬ 
pendence to the exponential time decay we used N r = 2, 
3, 4, and 5, corresponding to N = 9, 13, 17, and 21 sites. 

Again, we generated each a(L , t) from adjusting late- 
tirne profiles to Eq. (55) , by allowing both a and b there 
to vary. We saw that the fitted value of b generally stayed 
between 0.15 and 0.28. 

We then fitted sequences of varying-A data for a(L , t) 
to the n = 1 term of Eq. (l25l) . i.e., o(A, f) = 
ao exp(— cq(L) t), with co(A) = c + d/L z . 

Allowing z to vary freely gave a large amount of scatter 
(0.5 < z < 3.5) among fits of four-A data for the three 
different geometries [chains, and nanotubes with either 
N w = 14 or unit aspect ratio]. We then recalled that, for 
ID systems in the low-current phase a<l/2or / 3<l/2 
(except on the coexistence line a = /3 < 1/2) the effec¬ 
tive exponent governing the approach to steady state is 
z' = 1 [13j. This is in contrast to the result from a rigor¬ 
ous Bethe ansatz calculation |23], namely z = 0, and can 
be explained by a mean-field continuum formulation re¬ 
lated to kinematic-wave propagation jl3|. Thus we plot¬ 
ted our data for co(A) against 1/A, i.e. keeping z' = 1 
fixed. The results are shown in Fig. [8] It is seen that 
the numerical data for the sequences of co(A) fall reason¬ 
ably well onto a straight line consistent with z = 1, for all 
three geometries considered. From the vertical axis inter¬ 
cepts one gets respectively c = 0.009(5) (ID), c = 0.02(1) 
(nanotube with N w = 14), and c = 0.0046(44) (nanotube 
with unit aspect ratio). These are all definitely much 
lower than the mean-field prediction A 2 /2 = 0.08. It 
seems plausible from these data that the gap will vanish 
for very large nanotubes with finite aspect ratio (remain¬ 
ing finite in the quasi- and strictly ID cases). However, 
the relatively poor quality of the fits [ = 0.36, 10 and 

0.14, listed for each geometry in the same order as the c 
values ] indicates that a statement of this sort would have 
to be tested more extensively. 


C. p ^ 2q 

Initially we investigate the nanotube with p = q = 1 
at a c = 2(-\/2 — 1), /? c = 2 — \/2. These rates satisfy the 
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Figure 8. For systems with a = 0.3, /? = 0.7, plots 
of co(L) against 1/L, with co(L) defined via a(L,t) = 
ao exp(—co(L)t), the a(L,t) being given by fitting Eq. (l56l) 
to late-time profiles. Upper diagram: ID systems. Lower 
diagram: nanotubes with p = 2q = 1; squares: fixed width 
N w = 14; circles: fixed aspect ratio (A. R.)=l. 

conditions specified in Eqs. m urn for which the Mo- 
bius mapping predicts uniform steady state densities on 
each sublattice, though in general they remain distinct; 
namely, in this case they are x* = 2 — \/2, y* = — 1. 

For comparison, we consider also the chain with alter¬ 
nating bond rates p , 2q with p = q = 1/2. The mean-field 
Mobius mapping for this case coincides with that for the 
p = q = 1 nanotube, provided the injection/ejection rates 
are suitably renormalised, i.e. a = y/2— 1, /? = 1 — \/2/2. 
The respective steady state sublattice densities are then 
predicted to coincide, though of course the rate of ap¬ 
proach to steady state on the alternate-bond chain is half 
that for the nanotube. 

We took N w = 14, N r = 10 for the nanotube, and 
TV = 41 sites for the chain so both have the same number 
of sites along the flow direction. For the remainder of this 
Section, in both cases we always started with an empty 
lattice. 

Fig.[5]shows that the mean field prediction of flat sub¬ 
lattice density profiles in steady state is not fulfilled in 
numerical simulations. Also, the sublattice profiles for 
the nanotube and the alternating-bond chain do not co¬ 
incide, at variance with the fact that they share the same 
description via mean-field mapping. However, the mean 
field mapping predicts the steady-state sublattice den¬ 
sities to within at most 4% (for x *) or 16% (for y *) of 
numerical results. Since the predicted densities are them¬ 
selves separated by just over 40%, one can unequivocally 
ascribe each predicted sublattice profile to the correct 
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Figure 9. Steady state sublattice densities against position 
along flow direction for nanotube p = q = 1 (squares: x n t, 
y n t) at a = 2 ( y /2 — 1), (3 = 2 — %/2, and for chain with 
alternating bond rates p, 2 q with p = q = 1/2 (circles: x\ d, 
yid) at a = y/2 — 1, /3 = 1 — \/2/2 (see text). Horizontal 
dashed lines show mean field predictions applying for both 
cases: x* = 2 — \/2 , y* = \/2 — 1. 

numerically-generated subset of results. 

We defer further discussion of such discrepancies, and 
others which also pertain to steady-state aspects, to 
Sec. 1III Dl below. For the moment we investigate, for 
nanotubes with p = q = 1, the very late time behavior 
of the density profiles. Allowing for the observed non¬ 
uniformity of their limiting steady-state shapes, Eqs. prp 
and ( 1551 ) for a system at criticality should translate into: 

x t (t)=x}+G' sm — e~ x -^/ L)t (57) 

±J 

irf 

ye(t)=y* e +H' sin^e"*-^* (58) 

where now the position-dependent y/ are to be nu¬ 
merically obtained from steady-state simulation data. 

Results for the difference profiles, Sxe(t) = xe(t) — x\ 
and the similarly defined Sye(t), for the nanotube with 
p = q = 1, a = 2(^/2 — 1), ft = 2 — \/2 are exhibited in 
Fig. [TUI Late-time data were taken at t = 500 (for com¬ 
parison, the corresponding steady-state densities shown 
in Fig. [9] were taken at t = 2500). 

It is seen that the spatial dependence of 6xe(t) 
and 5ye(t) is indeed very close to that anticipated in 
Eqs. ( 1571 ) . ( 1551 ) . although the numerical results show a 
slight skew. The fit to a sine form shown as a dashed line 
in Fig. [Till corresponds to \| of = 49, which is unsatisfac¬ 
tory. We then allowed for a nonzero gap, by returning to 
the more general expressions Eqs. ( 1371 ) and ( 1551 ) . Fitting 
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to the n = 1 term of Eq. m, i-e., 


Sxe(t) 


—a{t) sin 



(59) 


we found the full-line curve depicted in Fig. 1101 with 
<j) = —0.022(1), x\ 0 i = 3-2. The small, but definitely non¬ 
zero, estimate of cj> is in line with the steady-state results 
shown in Fig. [3 in that both indicate the approximate, 
rather than exact, character of the mean-field description 
for p ^ 2 q. 

Furthermore, the difference profiles are almost entirely 
sublattice-independent, a feature which is not obviously 
forthcoming from the theory of Sec. Ill Bl It can be shown 
(see Appendix El that this results from the existence 
of two distinct relaxation rates: one which is very fast, 
size-independent [which brings the sublattice profiles to 
shapes rather close to their steady-state ones] and a 
slower one, with characteristic times of the usual L z form. 
In the (not very short)-time regime for which the latter 
applies the sublattice distinction disappears for differ¬ 
ence profiles, and the dynamics can be described in an 
effective continuum approximation through linear equa¬ 
tions resulting from a Cole-Hopf transformation. For ex¬ 
ample, difference profiles taken at t = 250 for the sys¬ 
tem considered in Fig. [TU] already exhibit a degree of 
sublattice-independence very similar to that shown in 
the Figure. Using Eqs. dSZD, USED for simplicity, defin¬ 
ing G"{t) = G' e~ x -^/ L ' ,t one finds by fitting numeri¬ 
cal data G"(250)/G ,, (500) « 4.4, which corresponds to 
\-(ir/L) pa 6 x 10~ 3 . Direct evaluation via the theoreti¬ 
cal prediction Eq. 63), using the mean field values for x *, 
y*, r from Eqs. (1751) and (1571) gives A_(7r/L) = 1.0 x 10” 3 . 

Turning now to non-critical systems, proceeding 
along the lines followed above one can again adapt 
Eqs. (1371) . (1381) to make allowance for the position de¬ 
pendence of steady state profiles, for systems away from 
criticality but with a and /3 obeying Eq. (E3). 

For a = 0.4, /3 = 0.8 the numerically-obtained steady 
state profiles turned out to be nearly flat down to 3 — 4 
parts in 1000, with xe ~ 0.324, ye ps 0.209, except very 
near the system’s ends. These values are rather close to 
the mean field ones predicted via Eq. (1751) . namely xe = 
1/3, ye = 1/5 . The late-time difference profiles obtained 
in the way described above, at t = 100, are displayed in 
Fig. [TT] Fitting to Eq. (l5iH) gives a fairly good account 
of the behavior of 5xe(t) against t. also, the sublattice 
independence of difference profiles is obeyed to a good 
extent, though some slight discrepancies remain near the 
ejection end. From Eqs. (ESI, m CEB- theory predicts 
that the coefficient <f> in the position dependence of late¬ 
time density profiles should be <j> = (ln8)/2 = 1.04..., 
and that for the time dependence the slowest decay rate 
should be Xj = 0.166_ 

The fitting curve shown in Fig. [TT] corresponds to 
(f> = 0.34(1). A measure of self-consistency of the lat¬ 
ter can be gained by pointing out that, if c/)L >5 — 6 the 
minimum of Eq. (l59ll is located at £ ps L — (1/</>). Vi¬ 
sual inspection of Fig. [TT] confirms that numerical data 
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Figure 10. Late-time difference profiles, Sxe{t) = xe(t) — x\, 
and similarly for 5ye(t), against position along flow direction 
for nanotube p = q = 1 at a c = 2(\/2 — 1), /3 C = 2 — %/2, 
for t = 500. The dashed line is the fit of Sxe(t) to a sine 
form, see Eqs. 63, EHD- The full line is a fit of 5xe(t) to a 
sine-plus-exponential form, see Eq. El and text. 


indeed behave in this way. On the other hand, the mis¬ 
match between predicted and observed values of <f> is a 
rather extreme illustration of the limitations of mean field 
mapping predictions for p ^ 2 q, already evident e.g. in 
the density profiles of Fig. [3 

We checked the theoretical prediction for Ai by com¬ 
paring difference profiles at t = 80 with those for t = 100. 
Referring to Eq. (l5i?l) . one gets a(100)/a(80) = 0.08 ± 
0.05, broadly compatible with e ~ 20X T = 0.03615 .... 


D. Factorization in steady state 

It was seen in Sec. wn that numerical results 
for steady state density profiles on nanotubes and 
alternating-bond chains with p ^ 2q are at variance 
with the predictions of mean field Mobius mapping. Mis¬ 
matches of similar order have been found between mean- 
field results and numerical work regarding steady-state 
currents in graphene-like structures with p ^ 2q [l| . 

In the following, we expand on comments made in 
Ref.[l], regarding the issue of factorization in steady state. 

It is known for the strictly one-dimensional TASEP 
that, along a + /? = 1 the correlations vanish, i.e., the 
probabilities for occupation variables on different sites 
factorize [5]. As a consequence of this, along that line 
the mean field mapping produces exact results. For nan¬ 
otubes one can then check for factorization (or its ab¬ 
sence) , in order to test the extent to which the predictions 
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Figure 11. Late-time difference profiles, 8xe(t) = xe(t) — x^, 
and similarly for Sye(t), against position along flow direction 
for nanotube p = g= lata = 0.4, /3 = 0.8, for t = 100. The 
full line is the fit of 5xe(t) to a sine plus exponential form, 
using only the n = 1 term of Eq. (E3 [ in an adapted form to 
allow for the position dependence of steady state profiles, see 
text]. 


given via Mobius mappings are accurate. 

A direct test can be implemented by considering the 
(connected) correlation function, 

Cij = (Jij) -pij in) (1 - (Tj)) , ( 60 ) 

where (Jij) is the average current across a chosen bond ij 
with rate pij, connecting sites i, j with respective mean 
occupations (rj, (tj). Factorization then corresponds to 
Cij = 0 for all bonds ij. 

We have found that for the nanotube with p = 1, 
q = 1/2 Cij vanishes to the accuracy of simulation (typ¬ 
ically 1 part in 10 5 ) on (and only on) the line a + /3 = 1, 
the same as in the strictly one-dimensional case. This is 
a non-trivial higher-dimensional generalization of a well 
known result for the linear chain. On the other hand, 
with p = 1 = q = 1, we followed the predicted factor¬ 
ization line, Eq. (1261) . and found that in simulations of 
similar accuracy, the factorization is no better than 1 
part in 10 2 . This is illustrated in Fig. [T2J where data 
taken at the respective predicted critical points, namely 
a = /3 = 1/2 [p = 2q = 1] and a = 2(^2-1), /3 = 2-y/2 
[p = q = 1 ] are shown. For p = q = 1, data are shown 
also for (a,/3) = (0.4,0.8), i.e., further along the pre¬ 
dicted factorization line Eq. (1261) . 

Still with p = q = 1 we thoroughly scanned the (a, (3) 
parameter space, and found no evidence either of uniform 
sublattice profiles or of vanishing of Cij . 



I 


Figure 12. Nanotube with N w = 14, N r = 10: Cy of eq. (16011 . 
averaged over transverse coordinate, against position along 
flow direction. Full symbols: x— sublattice. Empty symbols: 
y— sublattice. For p = q = 1, ( a c , /3 C ) = (2(%/2 — 1), 2 — y/2). 


IV. DISCUSSION AND CONCLUSIONS 

We have presented a mean-field theory for the dynam¬ 
ics of driven flow with exclusion in graphene-like struc¬ 
tures, and numerically checked its predictions. 

For the special combination of bond rates p = 2q in 
the nanotube geometry, Eqs. © © show that the sub¬ 
lattice distinction goes away in mean field. So a con¬ 
tinuum picture can apply, giving Eq. CD for which a 
time-dependent solution is found by using the Cole-Hopf 
transformation. 

For the special boundary rates a = j) = 1/2 which cor¬ 
responds to criticality in the ID chain with uniform rates, 
predictions for the early- and late-time behavior of den¬ 
sity profiles are made respectively in Eqs. (1221) and (ET1) . 
These are borne out by numerics to very good accuracy, 
see Figs. [2] and [3] We focused on late-time behavior, for 
both ID and nanotube geometries, and showed that by 
systematically analyzing the results of density profile fits 
to Eq. (HHl) it was possible [see Eqs. (ICTl) (©2)] to ex¬ 
tract rather accurate estimates of the dynamic exponent 
z. For strictly ID systems, we find z = 1.51(1), in excel¬ 
lent agreement with the anomalous value z = 3/2 which 
is known Hid to apply in that case. For nanotubes, 
we found strong indications (see Fig. © that the limiting 
behavior for very long length depends on whether one 
considers (quasi-ID) systems of fixed width, or square¬ 
like ones with constant aspect ratio; while the former ex¬ 
hibit z again close to 3/2, the latter are characterized by 
z consistent with the mean-field value of 2 (within error 
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bars). In the standard language of critical phenomena, 
this would mean that the upper critical dimensionality 
for TASEP dynamics is certainly D c < 2. 

On the factorization line a + /3 = 1 where steady- 
state profiles are uniform both for uniform-rate chains 
and nanotubes with p = 2q [lj, we took a = 0.3, away 
from criticality. The main distinguishing feature here, 
relative to the critical case, is the opening of a gap of 
amplitude A 2 /2 = (1 — 2a) 2 /2, associated with the char¬ 
acteristic length A -1 . The predicted effects of this on 
late-time profile shapes are spelt out in Eq. (li?5l) . which 
is qualitatively supported by numerical data (see Fig. □ . 

However, the quantitative effects, on the density pro¬ 
files, of having A ^ 0 are not accurately described by 
the present mean-held theory. Partly because of this, 
attempts to extract the dynamical exponent z, by proce¬ 
dures similar to those followed in the gapless case, met 
with the difficulties described in Sec. ME 

We then resorted to an overall consistency check, based 
on keeping fixed the effective exponent value z' = 1 which 
holds for the low-current phase in ID systems |l3). The 
resulting fits of numerical estimates of the coefficients ap¬ 
pearing in the exponential time decay factor of Eq. (ESJ, 
shown in Fig. [51 produce a reasonably self-consistent pic¬ 
ture. 

For nanotubes with p ^ 2q (and chains with alternat¬ 
ing bonds), Fig. [5] illustrates that predictions for steady 
state profiles from mean held mapping are not as accu¬ 
rate as for p = 2 q, or for uniform chains. In particular, 
numerically-generated prohles display a distinctive de¬ 
gree of nonuniformity along the predicted factorization 
line. 

Since dynamics concerns the evolution from initial to 
steady state, rather than the detailed (time-independent) 
properties of the latter, we adapted our original for¬ 
mulation to allow for the observed non-uniformity of 
the sublattice-dependent limiting prohle shapes, see 
Eqs. (1571) . (1551) for critical systems, and Eq. (1551) for the 
off-critical case. We found that for late times the differ¬ 
ence prohles thus defined behave in a very close way to 
that predicted by the theory of Sec. Ill Bl see respectively 
Figs. M and Ef 

This latter remark deserves to be qualified, inasmuch 
as it refers strictly to the functional forms displayed in 
Eqs. (1571) . (1551) . or Eq.([55]) [respectively sine, or sine plus 
exponential ] rather than to numerical values of the asso¬ 
ciated parameters [respectively G", H ', or a(t), </>] which 
we estimate via best-fitting procedures. Although this is 
not as stringent a test of mean held theory as would be 
the case if the theory-predicted parameter values were 
used, working this way allows one to separate possible 
shortcomings of the mean held approximation in func¬ 
tional forms versus those in parameters. Furthermore, 
one can have a quantitative estimate (through \ 2 values) 
of discrepancies in mean held functional forms, rather 
than the qualitative impressions from the comparisons 
with full predictions coming from theory; one can also 
get quantitative estimates of parameters affected by fluc¬ 


tuation effects absent from mean held theory, with the 
hope that modest generalizations (like domain wall the¬ 
ory) might more accurately provide such parameters. 

An additional feature of the late-time difference 
prohles is that they are almost entirely sublattice- 
independent. This property has been shown (see Ap¬ 
pendix \M to result from the coexistence of two distinct 
relaxation rates: a very fast, size-independent one, and a 
slower one with characteristic times of the usual L z form. 
The latter applies, within an effective continuum picture, 
to the Goldstone modes resulting from particle number 
conservation. If one accepts that an accurate description 
of TASEP via mean held mapping goes together with full 
applicability of a continuum approximation, this would 
then explain why the late-time density differences gener¬ 
ally fall in line with mean-held, continuum-like, predic¬ 
tions. 

Detailed comparison of theoretical predictions from 
Sec. Ill Bl to numerical results beyond overall prohle 
shapes turns out to not be as accurate as for p = 
2 q. For the system considered in Fig. [Till theory gives 
for the exponential time-decay coefficient of Eq. 
X-(n/L) = 1 x 10~ 3 , while adjusting to numerical data 
gives A_(7 t/L) « 6 x 10~ 3 . Similarly, for the non-critical 
system corresponding to Fig. [TTJ using the theoretical 
prediction for Ai of Eqs. m, <GHD would give a ratio of 
difference-profile coefficients at t = 100 and t = 80 equal 
to 0.03615..., while this same ratio is estimated from 
numerical data as 0.08 ± 0.05. 

Finally, in Sec. MIDI we showed that a direct test of 
factorization of correlation functions in steady state pro¬ 
duces a clear correspondence between uniformity of ob¬ 
served steady state prohles, on the one hand, and numer¬ 
ical evidence of vanishing of correlations, on the other. 
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Appendix A: Fast transient equalization of 
sublattices 

For our purposes here, it is convenient to adopt the 
following notation: sites on the x— sublattice have even 
site label, with mean occupation xit', for the y— sublat¬ 
tice, with odd labels one has the mean occupation j/ 2 t+i- 
Bond rates are p and p' = 2 q. 

Thus the mean held defining equations for currents and 
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occupations, and their evolution, Eqs. © ©, become 


JitiL+i = px 2 e (1 — V 2 t+l) (Al) 

K 2 Z- 12 Z = p' 2 / 21-1 (1 — %2 z) ■ (A2) 

X2i = K 2 e~i2i — J 2 Z 2 Z +1 (A3) 

2/21+1 = J2t 21+1 — K2t+12t+2 ■ (A4) 

Eqs. (f7\3il) and ([Mil give: 
d 

-Q^ {x2i + V2Z+l) = A 21-121 — A 21+121+2 ■ (A5) 


When a continuum picture applies, the right-hand side of 
Eq. (EU becomes like a space derivative of K and is then 
small, so X 2 e + 2 / 21+1 becomes a slow variable; similarly 
for y 2 z~i + X2e- 

On the other hand, any linear combination a X 2 Z + 
b 2 / 21+1 with a ^ b decays rapidly towards zero. This 
implies that the rapid decay is towards "adiabatic" val¬ 
ues of x 2 e, 2 / 21+1 such that all I< 2 £-\ 21 - T 2121+1 and 
J2i2t+\ - A' 2 i+ 121+2 are zero. That is, 


AT 2 1-121 — J2Z2Z+1 — <72121+1 — AT 2 1+121+2 — • • • — C(t) . 

(A6) 

The function C[t ) is the adiabatically evolving "con¬ 
served current" related to the particle conservation rep¬ 
resented by the set of equations Eq. (IA5I) for all t. Those 
equations determine the adiabatic evolution of the con¬ 
served densities. 

After the very fast transients have died out the profiles 
on the two sublattices still differ from their steady-state 
values X 2 I 1 2 / 21+1 by amounts Sx 2 e(t), £> 2 / 21+1 (t); as shown 
in the following, such differences are essentially the same 


for either sublattice, as their approach to zero is governed 
by a single continuum-like evolution equation. 

The fast time scales for the evolution of ax 2 e + 5 2 / 21+1 
with a 7 ^ b, coming from equations without nearly can¬ 
celling currents, and so without conserved or spatial 
derivative aspects, have rates set just by p and p', and 
not by wave vectors or system size L. So they are of 
order one, rather than a power of L or wavelength. 

In the subsequent evolution (after the initial transient 
regime) (i) we can interpolate the density variables be¬ 
tween the sites of their sublattice, making very little er¬ 
ror; and (ii) use the resulting "continuumization" of sites 
to find the conserved current differences in terms of spa¬ 
tial derivatives: e.g., 2/21 is the interpolation of the odd 
sublattice variables t/ 21 - 1 , 2 / 21 + 1 ; similarly for x 2 z+i- So, 

A2I-I2I — K2l+12t+2 = 

= p' [y 2 e-i(l - X 21 ) - y 2 f+i(l - X 2 t+ 2 )\ ~ 

~p' [fad 1 - X2£+l)\ , (A7) 

and similarly for differences of adjacent J's. 

Combining Eqs. (IA5I) and (IA7I) [ and their counterparts 
for yu-i + x 2 e and J 2 z- 22 i-i - J 2 Z 21 + 1 , respectively], 
omitting the subscripts and tilde signs, redefining t as 
an "average" coordinate shared by a pair of adjacent x— 
and y— subllattice sites, and defining p{l) = ^(xz + yz ), 
one gets: 


dp 

m 



p { 1 - p) 


1 dp 

2 ~dl 


(A8) 


This is now the form which the Cole-Hopf transformation 
linearizes. 
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